Effects of network topology, transmission delays, and refractoriness on the response of 
coupled excitable systems to a stochastic stimulus 
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We study the effects of network topology on the response of networks of coupled discrete excitable 
systems to an external stochastic stimulus. We extend recent results that characterize the response 
in terms of spectral properties of the adjacency matrix by allowing distributions in the transmission 
delays and in the number of refractory states, and by developing a nonperturbative approximation 
to the steady state network response. We confirm our theoretical results with numerical simulations. 
We find that the steady state response amplitude is inversely proportional to the duration of re- 
fractoriness, which reduces the maximum attainable dynamic range. We also find that transmission 
delays alter the time required to reach steady state. Importantly, neither delays nor refractoriness 
impact the general prediction that criticality and maximum dynamic range occur when the largest 
eigenvalue of the adjacency matrix is unity. 



Networks of coupled excitable systems describe 
many engineering, biological and social applica- 
tions. Recent studies of how such networks re- 
spond to an external, stochastic stimulus have 
provided insight on information processing in sen- 
sory neural networks^ In agreement with re- 
cent experiments^, these studies showed that the 
dynamic range of neural tissue is maximized in 
the critical regime, which is precisely balanced 
between growth and decay of propagating exci- 
tation. This regime was studied theoretically for 
directed random Erdos-Renyi networks in Ref*i, 
where it was found to be characterized by a net- 
work mean degree equal to one. However, other 
studies^il showed that this condition does not 
specify criticality for other network topologies. In 
this paper, extending recent results, we present 
a general framework for studying the effects of 
network topology on the response to a stochastic 
stimulus. With this framework, we derive a re- 
quirement for criticality and maximum dynamic 
range that holds for a wide variety of network 
topologies. Moreover, we show that this predic- 
tion holds when refractory states and transmis- 
sion time delays are included in the network dy- 
namics, although other aspects of the response do 
depend on these properties. 



I. INTRODUCTION 

Many applications involve networks of coupled ex- 
citable systems. Two prominent examples are the spread 
of information through neural networks and the spread 
of disease through human populations. The collective 
dynamics of such systems often defy naive expectations 
based on the dynamics of their individual components. 
For example, the collective response of a neural network 



can encode sensory stimuli which span more than 10 or- 
ders of magnitude in intensity, while the response of a 
single neuron typically encodes a much smaller range of 
stimulus intensities. Likewise, the collective properties 
of social contact networks determine when a disease be- 
comes an epidemic. 

Recently, a framework to study the response of a 
network of coupled excitable systems to a stochastic 
stimulus of varying strength has been proposed. The 
Kinouchi-Copelli modeli considers the response of a di- 
rected Erdos-Renyi random network of coupled discrete 
excitable systems to a stochastic external stimulus. A 
mean-field analysis of this model predicted^ that the 
maximum dynamic range (the range of stimuli over which 
the network's response varies signicantly) occurs in the 
critical regime where an excited neuron excites, on aver- 
age, one other neuron. This criterion can be stated as the 
mean out-degree of the network being one, = 1, 

where the out-degree of a node d°"* is defined as the 
expected number of nodes an excited node excites in 
the next time step (Refii refers to this quantity as the 
branching ratio). 

Subsequent studies explored this system on networks 
with power-law degree distributions and hypercubic lat- 
tice coupling, and with a varying number of loops^ii^"— , 
showing that the criterion for criticality based on the 
network mean degree does not hold for networks with a 
heterogeneous degree distribution. However, these stud- 
ies (except^) do not take into account features that are 
commonly found in real networks, such as, for example, 
community structure, correlations between in- and out- 
degree of a given node, or correlations between the degree 
of two nodes at the ends of a given cdge^. Furthermore, 
they do not consider the effect of transmission delays or 
a distribution in the number of refractory states. 

In a recent reporlj^, we presented an analysis of the 
Kinouchi-Copelli model that accounts for a complex net- 
work topology. We found that the general criterion for 
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criticality is that the largest eigenvalue of the network 
adjacency matrix is one, A = 1, rather than = 1. 

While this improved criterion successfully takes into ac- 
count various structural properties of networks, our anal- 
ysis did not address the effect of delays or multiple re- 
fractory states, and was based on perturbative approxi- 
mations to the network response. In this paper, we will 
extend the results of Refi^ by developing a nonperturba- 
tive analysis that accounts for distributions in the trans- 
mission delays and number of refractory states. 

This paper is organized as follows. In Section|lI]we de- 
scribe previous related work and the standard Kinouchi- 
Copelli model. In Section IIIII we present the model to 
be analyzed and derive a governing equation for its dy- 
namics. In Section IIVI we present our main theoretical 
results. In Section |V] we apply our results to estimate 
the dynamic range of excitable networks. In Section IVII 
we present numerical experiments to validate our results. 
We discuss our results in Sec. IVIII 

II. BACKGROUND 

In this Section we describe the Kinouchi-Copelli 
modeli and other relevant previous work. In order to 
focus on the effects of network topology, the dynamics of 
the excitable systems is taken to be as simple as possible. 
The model considers coupled excitable elements. Each 
element i can be in one of m -I- 1 states, Xi. The state 
Xi = is the resting state, Xi ~ 1 is the excited state, and 
there may be additional refractory states Xi = 2, 3, m. 
At discrete times t = 0, 1, ... the states of the elements x* 
are updated as follows: (i) If element i is in the resting 
state, x* = 0, it can be excited by another excited ele- 
ment j, x* = 1, with probability Aij, or independently 
by an external process with probability 77; (ii) the ele- 
ments that are excited or in a refractory state, a;* > 1, 
will deterministically make a transition to the next re- 
fractory state if one is available, or return to the resting 
state otherwise (i.e., x-'*'^ = x- + 1 if 1 < < m, and 
= if a;* = m). 

For a given value of the external stimulation probabil- 
ity 77, which is interpreted as the stimulus strength, the 
network response F is defined in Reft^ as 

F^{f)u (1) 

where {■)t denotes an average over time and /* is the 
fraction of excited nodes at time t. Of interest is the 
dependence of the response F{'r]) on the topology of the 
network encoded by the connection probabilities Aij . In 
particular, it is found that, depending on the network A, 
the network response can be of three typcsii^,: quiescent, 
in which the network activity is zero for vanishing stim- 
ulus, lim _F 0; active, in which there is self-sustained 

activity for vanishing stimulus, lim F > 0; and critical. 

7J->-0 

in which the response is still zero for vanishing stimulus 
but is characterized by sporadic long lasting avalanches of 



activity that cause a much slower decay in the response, 
compared with the quiescent case case, as the stimulus 
is decreased. Recent experiments^ suggest that cultured 
and acute cortical slices operate naturally in the critical 
regime. Therefore, the network properties that charac- 
terize this regime are of particular importance. 

In Refji, the response F was theoretically analyzed 
as a function of the external stimulation probability, 77, 
using a mean-field approximation in which connection 
strengths were considered uniform, Aij = {d)/N for all 
i,j. It was shown that the critical regime is achieved 
at the value (d) = 1, with the network being quies- 
cent (active) if (rf) < 1 ((d) > 1). For more gen- 
eral networks (i.e., Aij not constant), (d) is defined as 
the mean degree (d) = iE^,J^^J = (d"') = (rf™*), 
where d-" = J2j ^ij ^'^d = Aji are the in- and 
out-degrees of node i, respectively, and (•) is an aver- 
age over nodes. Such critical branching processes result 
in avalanches of excitation with power-law distributed 
sizes. Cascades of neural activity with power-law size and 
duration distributions have been observed in brain tis- 
sue cultures^i2r22, awake monkeysiSiH, and anesthetized 
ratsi^"— . While {d) = 1 successfully predicts the critical 
regime for Erdos-Renyi random networks^, it does not re- 
sult in criticality in networks with a more heterogeneous 
degree distributioii^'i. Perhaps more importantly, previ- 
ous theoretical analysesi*^ii are not able to take into ac- 
count features that are commonly found in real networks, 
such as, for example, community structure, correlations 
between in- and out-degree of a given node, or correla- 
tions between the degree of two nodes at the ends of a 
given edge^. We will generalize the mean-field criterion 
(d) = 1 to account for complex interaction topologies 
encoded in the matrix A as well as refractoriness and 
transmission delays. 



III. GENERALIZED KINOUCHI-COPELLI 
MODEL 

A. Description of the model 

We will analyze a generalized version of the Kinouchi- 
Copelli model which includes possibly heterogeneous dis- 
tributions of delays and refractory periods. The model is 
as follows: 

• There arc TV excitable elements, labeled i — 
l,...,iV. 

• At discrete times t = 0, 1, each element i can be 
in one of -|- 1 states, x\. The state x\ = Q is the 
resting state, x* = 1 is the excited state, and there 
may be additional refractory states = 2, 3, to^. 

• If clement i is in the resting state at time t, x\ = 0, 
it can be excited in the next time step, x*"^^ ~ 1, 
by another excited element j with delay Tij (i.e., if 
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1) with probability Aij, or independently where r* in Eq. ([2]) is the rate of transitions from the 



by an external stimulus with probability 77. 

The elements that are excited or in a refractory 
state, > 1, will deterministically make a tran- 
sition to the next refractory state if one is avail- 
able, or return to the resting state otherwise (i.e., 
a;*+i = X* -t- 1 if 1 < .T* < m,, and x*+^ = if 



m, 



ready to the excited state, given by 



E 



A ■ ■ T*^'^'^ 



(7) 



• The coupling network, encoded by the matrix with 
entries Aij, is allowed to have complex topology. 

B. Model dynamics 



where /* is one if node j is excited at time t and zero 
otherwise, and E[-] denotes an ensemble average. As- 
suming that the neighbors of node i being excited are 
independent events, we obtain, letting = p*, 



By considering a large ensemble of realizations of the 
above stochastic process on the same network, we can 
define the probability that node i is at state x* at time t 
as pI{x). The probabilities evolve in one time step by 



p*+i(2)=p*(l). 



and we also have the normalization condition 



pm 



(2) 
(3) 
(4) 
(5) 



(6) 



(8) 



We note that the assumption of independence is reason- 
able if there are few short loops in the network, and has 
been successfully used in similar situation a^^i^^ . How- 
ever, this assumption is violated if the number of bidi- 
rectional links is significant, and therefore we will restrict 
our attention to purely directed networks. Inserting the 
expression above in Eq. ^ and eliminating in terms 
of p* for j = 2, . . . , m^, we obtain the governing equation 
for the dynamics of p* 



k=l / 



r/ + (1 - 77) 



TV 



(9) 



In the following, we will analyze the response of the net- 
work by studying solutions of this equation as a function 
of 77. 



IV. ANALYSIS 

In this Section we study the solutions of Eq. (9) and 
the associated network response. In Sec. IV A, we 
develop a nonpcrturbative approximation to the steady 
state response of the network. In Sec. IV B we analyze 
limiting cases of the steady-state response that give us 
additional qualitative insight. In Sec. IV C, we study 
the effect of a distribution in the transmission time de- 
lays on the time scale of relaxation to the steady state 
solutions. We then discuss in Sec. IV D how our results 
relate to previous work. 



A. Steady-state response 

First, we will study steady-state solutions to Eq. ([9]). 
To find those, we set p\ = pi in Eq. ([9]), which becomes 



Pj = (1 - m^Pi) 77 + (1 - 77) 



\-^[\-A,,p,) 



(10) 

Proceeding as in Ref.— , by assuming AijPj is small, we 
replace YI^{1 - A^jPj) by exp(- J2j ^ijPj) to get 



Pi^ {1- m^pi) {rj + (1 - 7/) 



(11) 



The assumption that Aijpj is small is motivated as fol- 
lows. If the weights Aij are not very different from 
each other and each node has many incoming connec- 
tions (such as in neural networks, where the number of 
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synapses per neuron is estimated^ to be of the order 
of 10,000), then near the onset of self-sustained activity 
one should have ^ij ~ 1 (the mean-field prediction of 
Refii, which we refine here, states that the node average 
of Aij is one at criticality), implying Aij is small. The 
quantity AijPj is even smaller, especially for low levels of 
activity where pj is small. 

To proceed further we find convenient to define an al- 
ternative network response F as 



where (.t) s substitute s w F{d) / (u) into 

Eq. ([TC]) yielding 



where 



F={f)t, 



(12) 



(13) 



and /* = 1 if node j is excited at time t and otherwise. 

The variable /* can be interpreted as proportional to the 
number of excited nodes weighted by their out-degree 
J2i Aij. In terms of the probabilities pi, F is 



^out 



F ^ 



(14) 



and can be interpreted as the fraction of links that suc- 
cessfully transmit an excitation. This is analogous to the 
interpretation of F in Eq. ([l} as the fraction of excited 
nodes. In principle, the definitions of / and F preclude 
their use in comparing directly against commonly used 
measures of activity since knowledge of the matrix A is 
required to estimate them. However, we note that in 
all the numerical experiments discussed below, F and F 
were found to be nearly identical. To develop a nonper- 
turbative approximation to F , we solve Eq. (|lll) for pi in 

terms of X]j=i^ijPj- Multiplying the resulting expres- 
sion by Afei and summing over j, we obtain 



AT 

E 

1=1 



N 



AkiPi = ^ Aki 



1 - (1 - /7)e- 



1 



(1 — 7]je ^j=i 



3 Pi 



(15) 

Now, we use the fact that the largest eigenvalue of A, A, 
is typically much larger than the second eigenvalu o^^'^° , 
and thus Ap w su, where s is a scalar to be determined, 
and u is the right eigenvector of A corresponding to A. 
The validity of this approximation will be discussed in 
Section IVICI With this substitution, the previous equa- 
tion reduces to 



SUk 



N 

E 

i=l 



A 



ki 



1 - (1 - ,j)e- 



1 + nii — mi(l — ry)e 



(16) 



Noting that 



F = 



^lout 



Fuk{d) 



N 



Y^Ak.- 

A 1 -L 



1 _ (1 _,y)e-F'".<'i)/W 



which may now be summed over k, simplified, and solved 
for F , yielding the scalar equation 



^out 1 — (1 — fi'j£-F''^{d)/{u) 

{d) 1 + m - m(l - 77)e--^"<'')/<") 



(17) 



We note that in the notation above, the outer average (•) 
corresponds to a sum over the index i in Eq. (|16p . Given 
the adjacency matrix A, Eq. (|17p can be solved numeri- 
cally to obtain the response F as a function of rj. We call 
Eq. ()17p the "nonperturbative approximation" since its 
derivation does not rely on a perturbative truncation of 
the product term of Eq. (llOp . and we will numerically test 
its validity in Sec. IVIl where we will find that Eq. (flT)) 
can be a good approximation for all values of rj. In or- 
der to gain theoretical insight into how some features of 
the network topology and the distribution of the num- 
ber of refractory states affect the response, we will use 
Eq. P?|) to obtain analytical expressions for the response 
in various limits. 



B. Perturbative approximations 

While the nonperturbative approximation developed 
in the last section provides information for all ranges of 
stimulus, it is useful to consider perturbative approxima- 
tions, for example, to determine what is the transition 
point from quiescent to active behavior. We will obtain 
an approximation to F which is valid for small 77 and F. 
To do this, we expand the right hand side of Eq. (fT6|) 
to second order in s and first order in i] (as we will see, 
expanding to second order in s is necessary to treat the 
77 = case) obtaining 



SUk 
N 



(18) 









a- 




+ SUi 









Y^Ak^ 



Multiplying by the left eigenvector entry Vk and summing 
over k we obtain, using AkiVk ~ Xvi and rearranging, 



= i]X{v) - siiX{vu{l + 2m)) (19) 
+ (A - l)s{vu). 



In terms of F , this equation becomes 

F{d)i^\{vu{l + 2m)) (u) + (A - l)F{d){vu)(u) 



F^{d)' 



+ m] ) A = 'n\{v){u)' 



(20) 
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To find the transition from no activity, F = 0, to self- 
sustained activity, F > 0, for vanishing stimulus, we let 
77 — >■ 0^ in the previous equation to find 





(\-l){vu){u) 



A < 1, 
A > 1, 



(21) 



where the solution F ~ was chosen for A < 1 to satisfy 
_F > 0. This equation shows that the transition from a 
quiescent network to one with self-sustained activity has, 
if the response F is interpreted as an order parameter, 
the signatures of a second order (continuous) phase tran- 
sition. In addition, while the eigenvalue A determines 
when this transition occurs (at A = 1), its associated 
eigenvectors u and v determine the significance of the ob- 
served response past the transition. If (uw^) ^ {uv){u), 
for example, the response might be initially too small to 
be of importance. One aspect that was not considered in 
Refi^ is how the distribution of refractory periods affects 
the response. If the refractory periods to,; are strongly 
positively correlated with the product v.iuf, they can sig- 
nificantly increase the term {mvu'^) in the denominator, 
decreasing the response. This can be intuitively under- 
stood by noting that this amounts to preferentially in- 
creasing the refractory period of the nodes that are more 
likely to be active (as measured by the approximation 
Pi oc Ui valid close to the critical regime), thus removing 
them from the available nodes for longer times. 

The response F for small stimulus and response in 
Eq. ((2T|) agrees with the perturbative expression derived 
for F directly from Eq. ([9]) in Ref.— if nii = 1 and A — )> 1, 
and confirms the findings in Ref.— that the critical point 
is determined by A = 1. Henceforth, we will refer to net- 
works with A < 1 as quiescent, to networks with A > 1 
as active, and to networks with A = 1 as critical. 

The behavior of the system for high stimulation is also 
of interest. When rj = 1, node i cycles dcterministically 
through its nii + l available states, and so pt = {l+mi)~^. 
The question is how this behavior changes as rj decreases 
from 1. This information can be extracted directly from 
Eq. (jlip by linearizing around the solution 77 = 1, = 
(1 -(-TTii)^^. Setting J] = 1 — 5j], pi = pi~Spi with Si] ^ 1, 
6 Pi Pi, we obtain 



5pi 



(22) 



Thus, the response of the nodes to a decreased stimu- 
lus depends on a combination of their refractory period 
(which determines pi) and decays exponentially with the 
number of expected excitations from its neighbors. In 
terms of the aggregate response F, Eq. ((22|) becomes, 
after multiplting by A^i, summing over k and i, and nor- 
malizing. 



dF_ 
drj 



{d 



id) 



(23) 



C. Dynamics near the critical regime 

As in Ref.—, we will study the transition from no ac- 
tivity to self-sustained activity in the limit of vanishing 
stimulus by linearizing Eq. ^ around p* = for 77 = 0. 
Assuming p* is small, we obtain to first order 



JV 

Pi — 2^ ^i]Pj 



(24) 



Assuming exponential growth, p* = a^Wi, we obtain 



N 



(25) 



The critical regime, determined as the boundary between 
no activity and self-sustained activity as 77 — 0, i.e, be- 
tween the solution = being stable and unstable, can 
be found by setting a = 1, obtaining 



N 



(26) 



This implies that the onset of criticality occurs when 
A = 1 and in this case w = u. This conclusion agrees with 
the results in Refj^ and those in the previous section. Al- 
though the critical regime is not affected by the presence 
of delays or refractory states, the rate of growth (decay) 
a of perturbations for the active (quiescent) regime de- 
pends on the distribution of delays. To illustrate this, 
we consider the case when the network deviates slightly 
from the critical state, so that the largest eigenvalue of A 
is A = l + SX and has right eigenvector u, Au = {l + 5\)u. 
Expecting the solution w to Eq. (|25p to be close to u, we 
set Wi = Ui + Sui and a = 1 + /i, where the rate of growth 
fi is assumed to be small. Inserting these in Eq. (|25p . we 
get to first order 



fiu + Su ^ uSX + ASu — fiAu, 



(27) 



where the entries of the matrix A are given by Aij = 
AijTij. To eliminate the term Su, we left-multiply by the 
left eigenvector of A, v, satisfying v'^A — {1 + SX)!;"^ . 
Canceling small terms, we get 



5X 



v"^ Au 



(28) 



If the delay is constant, Tij — r, we obtain 



5X 



l + r' 



(29) 



and in this particular case a more general result can be 
obtained from Eq. which implies a = A^/^^"*""^^. 
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D. Relation to previous results 

Here wc will briefly discuss how our results for the crit- 
ical regime agree with previous work in particular cases. 
Correlations between degrees at the ends of a randomly 
chosen edge (assortative mixing by degree^) can be mea- 
sured by the correlation coefRcient 

p = (30) 

defined in Ref.— , with (•)e denoting an average over 
edges. The correlation coefficient p is greater than 1 
if the correlation between the in-dcgrcc and out-degree 
of nodes at the end of a randomly chosen edge is pos- 
itive, less than one if the correlation is negative, and 
one if there is no correlation. For a large class of net- 
works, the largest eigenvalue may be approximated^^ 
by A w In the absence of correlations, 

when p = 1, the largest eigenvalue can be approximated 
by A ~ If there are no correlations be- 

tween and d""* at a node [node degree correlations) 
or if the degree distribution is sufficiently homogeneous, 
then « (d)^ and the approximation reduces to 

A « {d). This is the situation that was considered in 
Refji, and thus they found that the critical regime was 
determined by (d) = 1. In the case of Refs4iii, with more 
heterogeneous degree distributions, A ~ {d™d°^^)/{d) ap- 
plied, which accounts for their observation that the crit- 
ical regime did not occur at (d) = 1. 

The situation encountered here is analogous to what 
occurs in the analysis of the transition to chaos in 
Boolean nctworksi^ and in the transition to synchroniza- 
tion in networks of coupled oscillators^!, where it is found 
that, instead of the mean degree, the largest eigenvalue 
is what determines the transition between different col- 
lective dynamical regimes. 

Other previous studies in random networks have also 
investigated spectral properties of A to gain insight on 
the stability of dynamics in neural networks^^ and have 
shown how A could be changed by modifying the distri- 
bution of synapse strengths^^. In addition, it has been 
shown recently that the largest eigenvalues in the spec- 
trum of the connectivity matrix may affect learning effi- 
ciency in recurrent chaotic neural networks^!. 



V. DYNAMIC RANGE 

We have studied the response of the network to stimuli 
of varying strengths. In particular, we studied in de- 
tail the response close to the critical regime. As has 
been previously noted!, this regime corresponds to the 
point where the dynamic range A is maximized. In our 
context, the dynamic range can be defined as the range 
of stimulus rj that results in significant changes in re- 
sponse F. Typically the dynamic range is given in deci- 
bels and measured using arbitrary thresholds just above 
the baseline (lim^^o F = Fq) and below the saturation 



(lim^^i F = Fi) values, respectively, as illustrated in 
Fig.[T]for the active network case Fq > 0. More precisely, 
the value of stimulus rjiow {ilhigh) corresponding to a low 
(high) threshold of activity Flow {Fhigh) are found and 
the dynamic range is calculated as 

A = lOlogiQ{riMgh/mow)- (31) 



A 




Fi 








Sigh 






Fbw 








Fo 


; Hiow 


Hhigh 



I ToiogJn) 
A ' 



FIG. 1: Schematic illustration of the definition of dynamic 
range in the active network case. The baseline and saturation 
values are -Fb and Fi, respectively. Two threshold values, 
denoted by Fiow and Fhigh, respectively, are used to determine 
the range of values of 77 defined as the dynamic range A. 

Using our approximations to the response _F as a func- 
tion of stimulus 77, we can study the effect of network 
topology on the dynamic range. The first approximation 
is based on the analysis of Sec. IIV Al Using Eq. (|17p . the 
values of 77 corrcponding to a given stimulus threshold 
can be found numerically and the dynamic range calcu- 
lated. 

Another approximation that gives theoretical insight 
into the effects of network topology and the distribu- 
tion of refractory states on the dynamic range can be 
developed as in Ref.—, by using the perturbative approx- 
imations developed in Sec. IIVBI In order to satisfy the 
restrictions under which those approximations were de- 
veloped, we will use Fhigh = Fi and Flow = Fq ^ 1. 
Taking the upper threshold to be Fhigh = is reason- 
able if the response decreases quickly from Fi , so that the 
effect of the network on the dynamic range is dependent 
mostly on its effect on Flow- Whether or not this is the 
case can be established numerically or theoretically from 
Eq. (j22p . and we find it is so in our numerical examples 
when mi are not large (see Fig. [5]). Taking rjhigh = 1 and 
Viow = V* we have 

A = -101ogio(77*). (32) 

The stimulus level 77 can be found in terms of F by solv- 
ing Eq. (|20p and keeping the leading order terms in F, 
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obtaining 

_ F^d)^ {vu^ + m)) - F{d)i\ - l){u){uv) 

This equation shows that as 77 — the response scales 
as F ^ 1] for the quiescent curves (A < 1), and as 
F 77^/^ for the critical curve (A = 1). We highlight that 
these scaling exponents for both the quiescent and critical 
regimes are precisely those derived in Rcfii for random 
networks, attesting to their robustness to the generaliza- 
tion of the criticality criterion to A = 1, the inclusion of 
time delays, and heterogeneous refractory periods . This 
is particularly important since these exponents could be 
measured experimentallji^. Using this approximation for 
rj* in p2p , we obtain an analytical expression for the dy- 
namic range valid when the lower threshold F* is small. 
Of particular theoretical interest is the maximum achiev- 
able dynamic range A.„iax for a given topology. It can be 
found by setting A = 1 in Eq. p3p and inserting the result 
in Eq. ([5^. obtaining 



A- 



max 



= An - 10 lot 



10 



(34) 



where Ao = — 20 log]^g(F*) > depends on the thresh- 
old F* but is independent of the network topology or 
the distribution of refractory states. The second term 
of Eq. (p4| suggests that a positive correlation between 
refractoriness m and eigenvector entries u and v will de- 
crease dynamic range, whereas a negative correlation will 
increase dynamic range. This prediction may be investi- 
gated in more depth in future publications. The second 
term of Eq. also suggests that an overall increase 

in the number of refractory states will lead to an overall 
decrease in dynamic range. This is in contrast with the 
result of Refi^, which found that there exists a tti > 
which maximizes dynamic range in two-dimensional ar- 
rays of neurons. We note that the assumption of inde- 
pendence used in deriving Eq. (|S]) is not valid for a two 
dimensional array. 



VI. NUMERICAL EXPERIMENTS 

We tested our theoretical results from section IIVI by 
comparing their predictions to direct simulation of our 
generalized Kinouchi-Copelli ixiodel described in section 
mil Simulation parameters were chosen specifically to 
test the validity of Eqs. (H?]), ^ and All 

simulations, except where indicated, were run with N ~ 
10^ nodes for T = 10'^ timcstcps, over a range of r/ from 
10-5 to 1. 



A. Construction of networks 

We created networks in three steps: first, we created 
binary directed networks, Aij G {0, 1}, with particular 



degree distributions as described below, forbidding bidi- 
rectional links and self-connections; second, we assigned 
a weight to each link, drawn from a uniform distribu- 
tion between and 1; third, we calculated A for the 
resulting network, and multiplied A by a constant to 
rcscalc A to the targeted eigenvalue^!. The two classes of 
topology considered for simulations were directed Erdos- 
Rcnyi random networks and directed scale-free networks 
with power law degree distributions, where we set the 
power law exponent to 7 = 2.5, and enforced a minimum 
degree of 10 and a maximum degree of 1000. Erdos- 
Renyi networksS^ were constructed by linking any pair 
of nodes with probability p = 15/iV, and scale- free net- 
works were constructed by first generating in-degree and 
out-degree sequences drawn from the power law distribu- 
tion described above, assigning those target degrees to N 
nodes, and then connecting them using the configuration 
model2^. In some cases an additional fourth step was 
used to change the assortativity coefficient p, defined in 
Eq. ([50]) . of a critical (i.e., with A = 1) scale-free network, 
making this network more assortative (disassortative) by 
choosing two links at random, and swapping their desti- 
nation connections only if the resulting swap would in- 
crease (decrease) p. This swapping allows for the degree 
of assortativity (and thereby. A) to be modified while 
preserving the network's degree distribution^iii^. 



B. Results of numerical experiments 

We first demonstrate the ability of the non- 
perturbative approximation to predict aggregate network 
behavior in a variety of conditions. Fig. [2] shows a multi- 
tude of simulations (symbols) with the predicted behav- 
ior of Eq. ([T7)) overlaid (lines). The cases considered in 
Fig. [2] include different combinations of topology, assor- 
tativity, largest eigenvalue A, delays, and number of re- 
fractory states. The number of refractory states rrii was 
chosen either constant, rrii = m, or randomly chosen with 
equal probability among {1,2,..., m,„ax}- Similarly, the 
delays were either constant, Tij = t, or uniformly 
chosen with uniform probability in (0,Tniax)- The pre- 
dictions capture the behavior of the simulations, with 
particularly good agreement for networks with neutral 
assortativity, p = I- In the assortative and disassortative 
cases shown [cases (3) and (6) in Fig. [5], low and high 
stimulus simulations are well captured by the prediction, 
while a small deviation can be observed for intermedi- 
ate values of rj [e.g., in case (6) in the right panel of 
Fig. [5J the crosses have a small systematic error around 
1] = 10^2] . In Sec. IVICi we wiU discuss why Eq. pT)) . 
which assumes Ap w sw, works so well. In particular, we 
will discuss why this approximation is expected to work 
well for small and large 77. 

As reported previously^, we find in our simulations 
that networks with A = 1 show critical dynamics and 
exhibit maximum dynamic range. This applies to ran- 
dom networks, scale free networks, and scale free net- 
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FIG. 2: Semi-log plots of data from five simulations (symbols) testing a variety of situations in order to show the robustness of 
Eq. (|17|l (lines) to various sets of conditions: (1) Random network; A = 1; mixed refractory states, rm £ {1, 2, 3}; no delays; (2) 
Random network; A = 0.7; no refractory states, rrii — m = 1; mixed delays, Tij £ {0, 1,2,3}; (3) Scale free network; A = 0.8; 
disassortative rewiring; no refractory states, rrii — m = 1; no delays; (4) Scale free network; A — 1.2; uniform refractory states, 
rrii = m = 2; no delays; (5) Scale free network; A = 1.0; mixed refractory states, rrii £ {1, 2, 3, 4}; mixed delays, Tij £ {0, 1, 2}; 
(6) Scale free network; A — 1.2; uniform refractory states, rrii = m — 1; no delays. The plots show excellent agreement between 
prediction and simulation at many points in parameter space. 



works with modified assortativity. Networks with A < 1 
exhibit no self-sustained activity in the absence of stim- 
ulus, whereas networks with A > 1 exhibit self-sustained 
activity. Furthermore, in all numerical experiments, with 
distributed refractory states, and various time delays, the 
criticality of networks at A = 1 was preserved as predicted 
above. Typical results in Fig. |3] (a) show the response F 
as a function of stimulus 77 for scale free networks with 
7 = 2.5, refractory states = m = 1, and no time de- 
lays, with A ranging from 0.2 to 1.8. Each symbol in the 
figure is generated by a single simulation on a single net- 
work realization. Lines show F obtained from numerical 
solution of Eq. PT)) . We note that the simulations with 
A = 1 show a deviation from the theoretically predicted 
critical curve for values of rj less than 10^^. We believe 
this is due to the fact that for such low values of 77, a much 
longer time average than the one we are doing would be 
required. For example, with 77 = 10~'^ we expect that, 
using 10^ time steps, a given node will not be excited 
externally with probability w 0.37. This might be 
especially important in the critical regime, where activity 
is mostly determined by sporadic avalanches propagated 
by hubs. 

Figure |3] (b) shows the dynamic range A calculated 
using F* = 10~^ directly from the simulation (circles) 
and using Eq. ((T7| (dashed line). As demonstrated in 
Refi^, the dynamic range is maximized when A = 1. We 
note that in Refi^ the dynamic range was estimated using 
a perturbative approximation, and as a consequence our 
prediction had a systematic error in the A > 1 regime [cf. 
Fig. 1 (b) in Refj^]. The nonperturbative approximation 
Eq. P7)) results in a much better prediction. 

Figure S] shows the transition that occurs at A = 1 
when 77 for experiments with a varying number of 
refractory states, m = 1, 3, and 5. Symbols indicate the 
results of direct simulation using 77 = 10~^, and the lines 
correspond to Eq. ([TT]), which describes well the result 



of the simulations. We found that for this particular 
network, the perturbative approximation (j2ip only gives 
correct results very close to the transition at A = 1, and 
its quantitative predictions degrade quickly as F grows. 
[A similar situation can be observed in Fig. 2(b) of Rcfi^.] 
However, we found that the perturbative approximation 
is still useful to predict the effect of the refractory states. 
Eq. (PT|) predicts that the response should scale as (m -I- 
1/2)^^. The inset shows how, after multiplication by 
(?Ti+l/2), the response curves collapse into a single curve. 
Figure |4] also depicts a linear relationship, F ~ (A — 
1) for A > 1. Making a connection with the theory of 
nonequilibrium phase transitions in which F ^ (A — Ac)'^, 
we derive Ac = 1 and the critical exponent /3 = 1. 

Figure [5] shows the response F close to 77 = 1 calculated 
for various values of tti from the simulation (symbols), 
and from Eq. (solid lines). Eq. describes well 
the slope of F close to 77 = 1. An important observation 
is that as m grows, the relative slope F~^dF /dr] at 77 = 1 
decreases. Therefore, if the typical refractory period m 
is large, the response F saturates [e.g., reaching 90% of 
F{1)] for smaller values of r/. 

Transmission delays, as in the analogous system of 
gene regulatory networksi^, do not modify steady state 
response. However, delays modify the time scale of re- 
laxation to steady state. We quantified this modifica- 
tion in the growth rate in Eq. (|28p . which determines 
the growth rate of perturbations from an almost critical 
quiescent network in terms of a matrix determined from 
the distribution of delays. In Fig. |6]we show time series 
(solid lines) for the initial growth in the number of ex- 
cited nodes within four active networks with and without 
time delays. For comparison, we show the slope that re- 
sults from the corresponding growth rates obtained from 
Eq. (|28p (dotted lines). The timesteps shown on the 
horizontal axis have been shifted to display multiple re- 
sults together, but not rescaled or distorted. As shown 
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FIG. 3; Simulation data for scale-free networks of 10* nodes 
(symbols) and numerical solution of Eq. 1171 (lines'), (a) Stim- 
ulus vs response predictions agree well in the regime where 
Ap « su, as discussed in section IVI CI Eigenvalues range 
from 0.2 to 0.9 (blue squares), exactly 1.0 (red diamonds), 
and from 1.1 to 1.8 (black circles), (b) Dynamic range predic- 
tions capture maximization at A = 1 as well as the non-critical 
behavior. 



in Fig. ini Eq. ()28|) is helpful in quantifying the growth 
rate of signals within the network in the regime during 
which growth is exponential. In this limited regime, sim- 
ulation data compare well with time series of excitations, 
and capture the growth rate's dependence on eigenvalue 
and time delays. We note here that Eq. ([25]) predicts the 
growth rate of p\ , and therefore the growth rate of both 
/* and /*. Here we have chosen to show the growth in 
the number of excited nodes (proportional to /*) which 
is more experimentally accessible than /*. 



C. Validity of the approximation Ap oc u 



0.15 




0.05 



largest eigenvalue, \ 



FIG. 4: Phase transitions of F,,->o for different refractory 
states, m for simulations (symbols) and Eq. (|17[l (lines). In- 
set: Eq. (|2ip predicts that phase transitions should scale by 
(m-f 1/2)^^, confirmed by rescaling data from the larger plot 
accordingly. 




FIG. 5: Simulation data (symbols) compare reasonably with 
the prediction of the perturbative approximation close to sat- 
uration, Eq. 1231 for different refractory states, drj was chosen 
to be the different between r) = lO" and rj = 10""'^, corre- 
sponding to the two rightmost data points of each simulation. 



Here we will address the question of the validity of 
our approximation Ap cx u, which was used to develop 
the nonperturbative approximation Eq. (|17p . First, we 
note that when and p are small, the linear analysis of 
Sec, live] and Refi^ shows that p (x u, and therefore the 
approximation Ap oc m is justified in this regime. As rj 
grows, and for situations where p is not small, one should 
expect deviations of p from being parallel to u. How- 
ever, we note that since pi measures how active node i 



is, it should still be highly correlated with the in-degree 
of node i. Since in many situations the in-degree is also 
correlated with the entries of the eigenvector u, we ex- 
pect that in those cases p remains correlated with u. Af- 
ter multiplication by A, the approximation can only be- 
come better. For the class of networks in which the ratio 
between the largest eigenvalue A and the next largest 
eigenvalue scales as -y/ (d) (which include Erdos-Rcnyi 
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FIG. 6: Time series (solid lines) for initial growth of signals 
within four active networks, with growth rates from Eq. (|23p 
shown (dotted lines). The timesteps shown on the horizontal 
axis have been shifted to display multiple results together, 
but not rescaled or distorted. In less than 100 timesteps, 
all networks tested exhausted the exponential growth regime. 
N=100000 nodes and r) = 10-^ for (1) A = 1.1, r = 0, (2) 
A = 1.2, r = 0, (3) A = 1.1, r G {0,1,2,3}, (4) A = 1.2, 
rG{0,l,2,3}. 



and other networks^S), we expect that Ap oc u should 
be a good approximation. 

Another reason why the approximation Ap oc u works 
well even when F is not smaU is that the errors in- 
troduced by this approximation vanish exactly when 
77 = 1. To see this, note that for 77 = 1, since each 
node cycles deterministically through its + 1 avail- 
able states, we have pi = 1/(1 + mi), which gives F = 
E,,A.,(l+m,)-7E,,^.j = (d°"*(l + m)-i)/(d), 
which agrees exactly with the result of setting 77 = 1 in 
Eq. p7p . Thus, even as the assumption Ap cx u may 
become less accurate as 77 grows, the importance of the 
error introduced by it decreases and eventually vanishes 
at 77 = 1. 

To illustrate how the assumption Ap « su works in 
some particular examples, Fig. [7] compares normalized 
Api and Ui for a variety of eigenvalues and stimulus lev- 
els. Good agreement between them (characterized by a 
high correlation) indicates that the assumption of section 
IIV Al is valid, whereas more noisy agreement for some 
cases indicates that the assumption Ap cx u is invalid 
(although, as discussed above, this does not necessar- 
ily imply that the nonperturbative approach will fail). 
Low stimulus levels in quiescent networks (top left panel) 
show relatively low correlation for short simulations, but 
the correlation improves with more timesteps as relative 
nodal response increases at well connected nodes and de- 
creases at poorly connected nodes. Assortative networks 
(bottom panels) show slightly lower correlation as well, 
corroborating the results shown in Fig. [5] where the pre- 
dictive power of Eq. (|17p is slightly diminished for the 
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FIG. 7: Plots of normalized Api vs sit, for scale-free networks, 
with eigenvalues 0.6 (blue, top row), 1.0 (red, middle row), 
and 1.2 with assortative mixing (black, bottom row) at stimu- 
lus levels Tj — 10"''', 10"'^, and 1 for the left, middle, and center 
columns respectively. Agreement is very good for critical and 
active cases, with more noise in the quiescent case due to less 
incoming stimuli over the duration of the simulation. 

assortative network. As expected, correlation between 
Ap and u entries is worst at 77 = 1 (right panels), but 
we reiterate that for rj ~ \ this error does not affect the 
predictions of Eq. p?]) . 



VII. DISCUSSION 

In this paper we studied a generalized version of the 
Kinouchi-Copelli model in complex networks. We devel- 
oped a nonperturbative treatment [Eq. ([T7| ] that allows 
us to find the response F of the network for a given value 
of the stimulus given a matrix of excitation transmission 
probabilities A. Our approach includes the possibility 
of heterogeneous distributions of excitation transmission 
delays and numbers of refractory states. An important 
assumption in our theory is that there are many incom- 
ing links to every node, which allows us to transform the 
product in Eq. ([9]) into an exponential. This assumption 
is very reasonable for neural networks, where the number 
of synapses per neuron is estimated^ to be of the order of 
10, 000. In addition, in order to obtain a closed equation 
for F, we assumed Ap oc u. As discussed in Sec. IVI CI 
this approximation works well in the regime when the 
response and stimulus are small. Furthermore, the er- 
ror introduced by this approximation becomes smaller as 
the probability of stimulus increases and eventually van- 
ishes for 77 = 1. The result is that Eq. ([T7)) predicts the 
response F satisfactorily for all values of 77. While we 
validated our predictions using scale-free networks with 
various correlation properties, we did not test them in 
topologies in which mean-field theories have been shown 
to fail, such as periodic hypercubes and branching tree 
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nctwork a^'^^ . This study is left to future research. 

Our theory describes how the introduction of addi- 
tional refractory states modifies the network response by 
modifying Eq. P7|) . In addition, their effect is captured 
by the perturbative approximations of Sec. IIVBI which, 
although valid in principle only for very small F, we have 
found successfully predict the effect of a distribution in 
the number of refractory states for a larger range of re- 
sponse values. 

We studied the effect of time delays on the time scale 
needed to reach a steady-state response, and found that 
Eq. determines the growth rate of perturbations 

from a quiescent, almost critical network. The tempo- 
ral characteristics of the response could be important in 
the study of sensory systems, in which the stimulus level 
might be constantly changing in time. Additionally, de- 
lays may be important in studying the phenomenon of 
synchronization and propagation of wavefronts, which we 
do not study here. Synchronization in epidemic mod- 
els similar to the model considered here has been well- 
described in the absence of time delays^^, and synchro- 
nization in Rulkov neurons has been shown to be affected 
subtly by time delays^. However, the effect of time de- 
lays on synchronization in our model remains an open 
line of inquiry. 

An important practical question regarding the applica- 
tion of our theory to neuroscience is how our results can 
be made compatible with the presence of excitatory and 
inhibitory connections in neural networks. Considering 
one excited neuron, and after excitatory and inhibitory 
connections are taken into account, the important quan- 
tity that determines the future activity of the network 



is how many other neurons are expected to be excited 
by the originally excited neuron. This number might de- 
pend on the overall balance of excitatory and inhibitory 
connections, but it must be a positive number. The 
Kinouchi-Copelli model we are using, and similar models 
used successfully by neuroscientists to model neuronal 
avalanches^, have therefore considered only excitatory 
neurons, while adjusting the probabilities of excitation 
transmission to account for different balances of excita- 
tory and inhibitory neurons. Nevertheless, we believe 
a generalization of the Kinouchi-Copelli model that ac- 
counts for inhibitory connections should be investigated 
in the future. 

Another important issue is the generality of our find- 
ings for more biologically realistic excitable systems. We 
conjecture that the effect of network topology on the dy- 
namic range of networks of continuous-time, continuous- 
state coupled excitable systems such as coupled ODE 
neuron models^ is qualitatively similar to its effect on 
the class of discrete-time, discrete-state dynamical sys- 
tems studied here. However, this remains open to inves- 
tigation. 
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